#############################################################
### Replication Code: Appendices B, C, D, E, and G        ###
### Title: Can Conservatives Be Persuaded to Support UBI? ###
### Author: Eddy S. F. Yeung                              ###
### Version: September 6, 2022                            ###
#############################################################

# Set-up ----
## Clean the R environment and set the working directory
rm(list = ls())
setwd("~/Desktop/UBI/PB_replication") # change to your own directory

## Load the required packages
library(tidyverse) # version 1.3.1
library(multcomp)  # version 1.4-17
library(broom)     # version 1.0.0
library(estimatr)  # version 0.30.4
library(texreg)    # version 1.37.5
library(gridExtra) # version 2.3
library(extrafont) # version 0.17
library(interflex) # version 1.2.6

## Import the cleaned dataset
load("cleaned.RData")

# Analysis for Appendix B ----
## Figure A4: Balance on Demographic Variables across Experimental Groups ----
### Create an empty PDF file to store the graph
pdf("fig_balance.pdf", height = 11, width = 11)
par(mfrow = c(3, 3))

### Age
plot(density(df[df$group == "0", ]$age, na.rm = T, bw = 4), 
     xlab = "Age", ylim = c(0, .025), xlim = c(10, 100), 
     main = "Age", lty = 1, cex = 1.2, family = "Times")
par(new = T)
plot(density(df[df$group == "1", ]$age, na.rm = T, bw = 4), 
     xlab = "", ylim = c(0, .025), xlim = c(10, 100), 
     main = "", axes = F, ylab = "", lty = 2)
par(new = T)
plot(density(df[df$group == "2", ]$age, na.rm = T, bw = 4), 
     xlab = "", ylim = c(0, .025), xlim = c(10, 100), 
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topright", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")

### Sex
plot(density(df[df$group == "0", ]$female, na.rm = T, bw = 0.1), 
     xlab = "Sex (0 = Male; 1 = Female)", ylim = c(0, 2.7), xlim = c(-0.4, 1.4), 
     main = "Sex", lty = 1, cex = 1.2)
par(new = T)
plot(density(df[df$group == "1", ]$female, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 2.7), xlim = c(-0.4, 1.4), 
     main = "", axes = F, ylab = "", lty = 2)
par(new = T)
plot(density(df[df$group == "2", ]$female, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 2.7), xlim = c(-0.4, 1.4),
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topright", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")

### Marital status
plot(density(df[df$group == "0", ]$married, na.rm = T, bw = 0.1), 
     xlab = "Married (0 = No; 1 = Yes)", ylim = c(0, 2.5), xlim = c(-0.4, 1.4), 
     main = "Marital Status", lty = 1, cex = 1.2)
par(new = T)
plot(density(df[df$group == "1", ]$married, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 2.5), xlim = c(-0.4, 1.4), 
     main = "", axes = F, ylab = "", lty = 2)
par(new = T)
plot(density(df[df$group == "2", ]$married, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 2.5), xlim = c(-0.4, 1.4),
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topright", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")

### Education
df$college <- ifelse(df$educ >= 4 & df$educ <= 7, 1,
                     ifelse(df$educ >= 0 & df$educ <= 3, 0, NA))
plot(density(df[df$group == "0", ]$college, na.rm = T, bw = 0.1), 
     xlab = "College Degree (0 = No; 1 = Yes)", ylim = c(0, 3), xlim = c(-0.4, 1.4), 
     main = "Education", lty = 1, cex = 1.2)
par(new = T)
plot(density(df[df$group == "1", ]$college, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 3), xlim = c(-0.4, 1.4), 
     main = "", axes = F, ylab = "", lty = 2)
par(new = T)
plot(density(df[df$group == "2", ]$college, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 3), xlim = c(-0.4, 1.4),
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topright", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")

### Income
df$income3 <- ifelse(df$income >= 0 & df$income <= 4, 0, 
                     ifelse(df$income >= 5 & df$income <= 7, 1, 
                            ifelse(df$income >= 8 & df$income <= 12, 2, NA)))
plot(density(df[df$group == "0", ]$income3, na.rm = T, bw = 0.25), 
     xlab = "Household Income (1 = At Medium Income; 0 = Below; 2 = Above)", 
     ylim = c(0, .85), xlim = c(-.7, 2.7), xaxt = "n", 
     main = "Household Income", lty = 1, cex = 1.2)
axis(1, at = c(0, 1, 2))
par(new = T)
plot(density(df[df$group == "1", ]$income3, na.rm = T, bw = 0.25), 
     xlab = "", ylim = c(0, .85), xlim = c(-.7, 2.7), 
     main = "", axes = F, ylab = "", lty = 3)
par(new = T)
plot(density(df[df$group == "2", ]$income3, na.rm = T, bw = 0.25), 
     xlab = "", ylim = c(0, .85), xlim = c(-.7, 2.7), 
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topright", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")

### Race
plot(density(df[df$group == "0", ]$white, na.rm = T, bw = 0.1), 
     xlab = "White (0 = No; 1 = Yes)", ylim = c(0, 3), xlim = c(-0.4, 1.4), 
     main = "Race", lty = 1, cex = 1.2)
par(new = T)
plot(density(df[df$group == "1", ]$white, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 3), xlim = c(-0.4, 1.4), 
     main = "", axes = F, ylab = "", lty = 2)
par(new = T)
plot(density(df[df$group == "2", ]$white, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 3), xlim = c(-0.4, 1.4),
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topleft", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")

### Employment status
df$employ2 <- as.numeric(df$employ) - 1
plot(density(df[df$group == "0", ]$employ2, na.rm = T, bw = 0.25), 
     xlab = "Employment (0 = Employed; 1 = Not Employed; 2 = Retired)", 
     ylim = c(0, 1), xlim = c(-.7, 2.7), xaxt = "n", 
     main = "Employment Status", lty = 1, cex = 1.2)
axis(1, at = c(0, 1, 2))
par(new = T)
plot(density(df[df$group == "1", ]$employ2, na.rm = T, bw = 0.25), 
     xlab = "", ylim = c(0, 1), xlim = c(-.7, 2.7), 
     main = "", axes = F, ylab = "", lty = 3)
par(new = T)
plot(density(df[df$group == "2", ]$employ2, na.rm = T, bw = 0.25), 
     xlab = "", ylim = c(0, 1), xlim = c(-.7, 2.7), 
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topright", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")

### Union membership
plot(density(df[df$group == "0", ]$union, na.rm = T, bw = 0.1), 
     xlab = "Household Union Membership (0 = No; 1 = Yes)", 
     ylim = c(0, 3.5), xlim = c(-0.4, 1.4), 
     main = "Union Membership", lty = 1, cex = 1.2)
par(new = T)
plot(density(df[df$group == "1", ]$union, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 3.5), xlim = c(-0.4, 1.4), 
     main = "", axes = F, ylab = "", lty = 2)
par(new = T)
plot(density(df[df$group == "2", ]$union, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 3.5), xlim = c(-0.4, 1.4),
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topright", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")

### Party identification
df$republican <- ifelse(df$pid >= 4 & df$pid <= 6, 1, 
                        ifelse(df$pid >= 0 & df$pid <= 3, 0, NA))
plot(density(df[df$group == "0", ]$republican, na.rm = T, bw = 0.1), 
     xlab = "Republican (0 = No; 1 = Yes)", ylim = c(0, 2.5), xlim = c(-0.4, 1.4), 
     main = "Party Identification", lty = 1, cex = 1.2)
par(new = T)
plot(density(df[df$group == "1", ]$republican, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 2.5), xlim = c(-0.4, 1.4), 
     main = "", axes = F, ylab = "", lty = 2)
par(new = T)
plot(density(df[df$group == "2", ]$republican, na.rm = T, bw = 0.1), 
     xlab = "", ylim = c(0, 2.5), xlim = c(-0.4, 1.4),
     main = "", axes = F, ylab = "", lty = 3)
par(new = T, family = "Times")
legend("topright", legend = c("Control Group", "Treatment Group 1", "Treatment Group 2"), 
       lty = c(1, 2, 3), cex = 1, bty = "n")
dev.off()

# Analysis for Appendix C ----
## Table A1: Coding of Conservative Answers of Policy Statement Questions (column 3) ----
mean(df$policy1, na.rm = T)
mean(df$policy2, na.rm = T)
mean(df$policy3, na.rm = T)
mean(df$policy4, na.rm = T)
mean(df$policy5, na.rm = T)
mean(df$policy6, na.rm = T)
mean(df$policy7, na.rm = T)
mean(df$policy8, na.rm = T)
mean(df$policy9, na.rm = T)
mean(df$policy10, na.rm = T)

## Figure A5: Distribution of Conservatism Scores ----
fig_cons.score <- 
  ggplot(data = df, aes(x = cons.score)) +
  geom_histogram(bins = 11, color = "black", fill = "grey") +
  scale_x_continuous(breaks = seq(0, 10, by = 2)) +
  xlab("Conservatism Score (0 = Least Conservative; 10 = Most Conservative)") +
  ylab("Number of Respondents") +
  coord_cartesian(xlim = c(-0.5, 10.5), ylim = c(0, 500)) +
  theme_bw() +
  theme(text = element_text(family = "Times", size = 13))
ggsave(file = "fig_cons.score.pdf", fig_cons.score, width = 8, height = 4)

# Analysis for Appendix D ----
## Percentage of self-reported conservatives in the sample ----
df$self.con <- ifelse(df$ideo == "4" | df$ideo == "5" | df$ideo == "6", 1, 0)
mean(df$self.con)

## Table A2: OLS Regression Corresponding to Equation 1 Based on a Self-Reported, Pretreatment Measure of Conservatism ----
### Specification 1
mod1.sicon <- lm_robust(support ~ treatment1 + treatment2 + self.con +
                          treatment1 * self.con + treatment2 * self.con,
                        data = df)

### Specification 2
mod2.sicon <- lm_robust(support ~ treatment1 + treatment2 + self.con +
                          treatment1 * self.con + treatment2 * self.con +
                          age + female + married + educ,
                        data = df)

### Specification 3
mod3.sicon <- lm_robust(support ~ treatment1 + treatment2 + self.con +
                          treatment1 * self.con + treatment2 * self.con +
                          age + female + married + educ + income + white + black +
                          employ + union + pid + welfare.receive,
                        data = df)

### Specification 4
mod4.sicon <- lm_robust(support ~ treatment1 + treatment2 + self.con +
                          treatment1 * self.con + treatment2 * self.con +
                          age + female + married + educ + income + white + black +
                          employ + union + pid + welfare.receive + welfare.bad +
                          racial.resent + demo.change,
                        data = df)

### Export the regression table to LaTeX
texreg(list(mod1.sicon, mod2.sicon, mod3.sicon, mod4.sicon),
       stars = c(0.01, 0.05, 0.10),
       include.ci = F,
       custom.header = list("Dependent Variable: Seven-Point UBI Support" = 1:4),
       custom.coef.names = 
         c("Constant", "Equalizing-Opportunity (EO) Frame",
           "Limiting-Government (LG) Frame", "Self-Reported (SR) Conservative (= 1)",
           "EO Frame × SR Conservative", "LG Frame × SR Conservative",
           "Age", "Female (= 1)", "Married (= 1)", "Education (7 = Highest)", 
           "Household Income (12 = Highest)", "White (= 1)", "Black (= 1)", 
           "Not Employed (= 1)", "Retired (= 1)", "Household Union Membership (= 1)", 
           "Party ID (6 = Strong Republican)", "Welfare Recipient (= 1)", 
           "Welfare System Evaluation (1 = Bad)", "Racial Resentment (8 = Highest)", 
           "Demographic Change (1 = Yes)"),
       custom.note = "Entries are OLS estimates with robust standard errors in parentheses.
       All significance tests are two-tailed with the following notations:
       $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$.",
       caption = "OLS Regression Corresponding to Equation 1 Based on a Self-Reported, Pretreatment Measure of Conservatism",
       fontsize = "small")

## Figure A6: Conditional Treatment Effects on Self-Reported Conservatives ----
### Create an empty data frame to store the estimates for later visualization
df.CATE1 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE1) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE1[1:3, 1] <- 1
df.CATE1[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (7-point scale DV)
### All respondents (including inattentive respondents)
H1a.full <- glht(mod3.sicon, linfct = c("treatment1 + treatment1:self.con = 0"))
df.CATE1[1, 2] <- confint(H1a.full)$confint[ , "Estimate"]
df.CATE1[1, 3] <- confint(H1a.full)$confint[ , "lwr"]
df.CATE1[1, 4] <- confint(H1a.full)$confint[ , "upr"]
df.CATE1[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3.alt1 <- lm_robust(support ~ treatment1 + treatment2 + self.con +
                         treatment1 * self.con + treatment2 * self.con +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                       data = subset(df, attention == 1 | attention == 2))
H1a.alt1 <- glht(mod3.alt1, linfct = c("treatment1 + treatment1:self.con = 0"))
df.CATE1[2, 2] <- confint(H1a.alt1)$confint[ , "Estimate"]
df.CATE1[2, 3] <- confint(H1a.alt1)$confint[ , "lwr"]
df.CATE1[2, 4] <- confint(H1a.alt1)$confint[ , "upr"]
df.CATE1[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3.alt2 <- lm_robust(support ~ treatment1 + treatment2 + self.con +
                         treatment1 * self.con + treatment2 * self.con +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                       data = subset(df, attention == 2))
H1a.alt2 <- glht(mod3.alt2, linfct = c("treatment1 + treatment1:self.con = 0"))
df.CATE1[3, 2] <- confint(H1a.alt2)$confint[ , "Estimate"]
df.CATE1[3, 3] <- confint(H1a.alt2)$confint[ , "lwr"]
df.CATE1[3, 4] <- confint(H1a.alt2)$confint[ , "upr"]
df.CATE1[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (7-point scale DV)
### All respondents (including inattentive respondents)
H2a.full <- glht(mod3.sicon, linfct = c("treatment2 + treatment2:self.con = 0"))
df.CATE1[4, 2] <- confint(H2a.full)$confint[ , "Estimate"]
df.CATE1[4, 3] <- confint(H2a.full)$confint[ , "lwr"]
df.CATE1[4, 4] <- confint(H2a.full)$confint[ , "upr"]
df.CATE1[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.alt1 <- glht(mod3.alt1, linfct = c("treatment2 + treatment2:self.con = 0"))
df.CATE1[5, 2] <- confint(H2a.alt1)$confint[ , "Estimate"]
df.CATE1[5, 3] <- confint(H2a.alt1)$confint[ , "lwr"]
df.CATE1[5, 4] <- confint(H2a.alt1)$confint[ , "upr"]
df.CATE1[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.alt2 <- glht(mod3.alt2, linfct = c("treatment2 + treatment2:self.con = 0"))
df.CATE1[6, 2] <- confint(H2a.alt2)$confint[ , "Estimate"]
df.CATE1[6, 3] <- confint(H2a.alt2)$confint[ , "lwr"]
df.CATE1[6, 4] <- confint(H2a.alt2)$confint[ , "upr"]
df.CATE1[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE1$group <- factor(df.CATE1$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE_self1 <- 
  ggplot(data = df.CATE1, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Self-Identified Conservatives") +
  ggtitle("Panel A: Seven-Point UBI Support") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.title = element_blank(),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-1.15, 1.15)) +
  guides(color = guide_legend(nrow = 1, byrow = T))
fig_CATE_self1

### Create an empty data frame to store the estimates for later visualization
df.CATE2 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE2) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE2[1:3, 1] <- 1
df.CATE2[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (binary DV)
### All respondents (including inattentive respondents)
mod3.binary0 <- lm_robust(percent.sup ~ treatment1 + treatment2 + self.con +
                            treatment1 * self.con + treatment2 * self.con +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = df)
H1a.binary0 <- glht(mod3.binary0, linfct = c("treatment1 + treatment1:self.con = 0"))
df.CATE2[1, 2] <- confint(H1a.binary0)$confint[ , "Estimate"]
df.CATE2[1, 3] <- confint(H1a.binary0)$confint[ , "lwr"]
df.CATE2[1, 4] <- confint(H1a.binary0)$confint[ , "upr"]
df.CATE2[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3.binary1 <- lm_robust(percent.sup ~ treatment1 + treatment2 + self.con +
                            treatment1 * self.con + treatment2 * self.con +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = subset(df, attention == 1 | attention == 2))
H1a.binary1 <- glht(mod3.binary1, linfct = c("treatment1 + treatment1:self.con = 0"))
df.CATE2[2, 2] <- confint(H1a.binary1)$confint[ , "Estimate"]
df.CATE2[2, 3] <- confint(H1a.binary1)$confint[ , "lwr"]
df.CATE2[2, 4] <- confint(H1a.binary1)$confint[ , "upr"]
df.CATE2[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3.binary2 <- lm_robust(percent.sup ~ treatment1 + treatment2 + self.con +
                            treatment1 * self.con + treatment2 * self.con +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = subset(df, attention == 2))
H1a.binary2 <- glht(mod3.binary2, linfct = c("treatment1 + treatment1:self.con = 0"))
df.CATE2[3, 2] <- confint(H1a.binary2)$confint[ , "Estimate"]
df.CATE2[3, 3] <- confint(H1a.binary2)$confint[ , "lwr"]
df.CATE2[3, 4] <- confint(H1a.binary2)$confint[ , "upr"]
df.CATE2[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (binary DV)
### All respondents (including unattentive respondents)
H2a.binary0 <- glht(mod3.binary0, linfct = c("treatment2 + treatment2:self.con = 0"))
df.CATE2[4, 2] <- confint(H2a.binary0)$confint[ , "Estimate"]
df.CATE2[4, 3] <- confint(H2a.binary0)$confint[ , "lwr"]
df.CATE2[4, 4] <- confint(H2a.binary0)$confint[ , "upr"]
df.CATE2[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.binary1 <- glht(mod3.binary1, linfct = c("treatment2 + treatment2:self.con = 0"))
df.CATE2[5, 2] <- confint(H2a.binary1)$confint[ , "Estimate"]
df.CATE2[5, 3] <- confint(H2a.binary1)$confint[ , "lwr"]
df.CATE2[5, 4] <- confint(H2a.binary1)$confint[ , "upr"]
df.CATE2[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.binary2 <- glht(mod3.binary2, linfct = c("treatment2 + treatment2:self.con = 0"))
df.CATE2[6, 2] <- confint(H2a.binary2)$confint[ , "Estimate"]
df.CATE2[6, 3] <- confint(H2a.binary2)$confint[ , "lwr"]
df.CATE2[6, 4] <- confint(H2a.binary2)$confint[ , "upr"]
df.CATE2[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE2$group <- factor(df.CATE2$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE_self2 <- 
  ggplot(data = df.CATE2, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr),
                position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Self-Identified Conservatives") +
  ggtitle("Panel B: Percentage Support for UBI") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-.3, .3)) +
  guides(colour = guide_legend(nrow = 1, byrow = TRUE))
fig_CATE_self2

### Combine the two plots (fig_CATE_self1 and fig_CATE_self2) and let them share a common legend
get_legend <- function(myggplot) {
  temp <- ggplot_gtable(ggplot_build(myggplot))
  leg <- which(sapply(temp$grobs, function(x) x$name) == "guide-box")
  legend <- temp$grobs[[leg]]
  return(legend)
}
legend <- get_legend(fig_CATE_self1) # save the legend
fig_CATE_self <- 
  grid.arrange(arrangeGrob(fig_CATE_self1 + theme(legend.position = "none"),
                           fig_CATE_self2 + theme(legend.position = "none"),
                           nrow = 1),
               legend, nrow = 2, heights = c(10, 1))
ggsave(file = "fig_CATE_self.pdf", fig_CATE_self, width = 8, height = 4.5)

# Analysis for Appendix E ----
## Percentage of moderates in the sample ----
df$moderate <- ifelse(df$cons.score == 3 | df$cons.score == 4, 1, 0)
table(df$moderate)
521 / (2009 + 521)

## Table A3: OLS Regression Corresponding to Equation 1 for Moderates ----
### Specification 1
mod1.moder <- lm_robust(support ~ treatment1 + treatment2 + moderate + 
                          treatment1 * moderate + treatment2 * moderate,
                        data = df)

### Specification 2
mod2.moder <- lm_robust(support ~ treatment1 + treatment2 + moderate +
                          treatment1 * moderate + treatment2 * moderate +
                          age + female + married + educ,
                        data = df)

### Specification 3
mod3.moder <- lm_robust(support ~ treatment1 + treatment2 + moderate +
                          treatment1 * moderate + treatment2 * moderate +
                          age + female + married + educ + income + white + black +
                          employ + union + pid + welfare.receive,
                        data = df)

### Specification 4
mod4.moder <- lm_robust(support ~ treatment1 + treatment2 + moderate +
                          treatment1 * moderate + treatment2 * moderate +
                          age + female + married + educ + income + white + black +
                          employ + union + pid + welfare.receive + welfare.bad +
                          racial.resent + demo.change,
                        data = df)

### Export the regression table to LaTeX
texreg(list(mod1.moder, mod2.moder, mod3.moder, mod4.moder),
       stars = c(0.01, 0.05, 0.10),
       include.ci = F,
       custom.header = list("Dependent Variable: Seven-Point UBI Support" = 1:4),
       custom.coef.names = 
         c("Constant", "Equalizing-Opportunity (EO) Frame",
           "Limiting-Government (LG) Frame", "Moderate (= 1)",
           "EO Frame × Moderate", "LG Frame × Moderate",
           "Age", "Female (= 1)", "Married (= 1)", "Education (7 = Highest)", 
           "Household Income (12 = Highest)", "White (= 1)", "Black (= 1)", 
           "Not Employed (= 1)", "Retired (= 1)", "Household Union Membership (= 1)", 
           "Party ID (6 = Strong Republican)", "Welfare Recipient (= 1)", 
           "Welfare System Evaluation (1 = Bad)", "Racial Resentment (8 = Highest)", 
           "Demographic Change (1 = Yes)"),
       custom.note = "Entries are OLS estimates with robust standard errors in parentheses.
       All significance tests are two-tailed with the following notations:
       $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$.",
       caption = "OLS Regression Corresponding to Equation 1 for Moderates",
       fontsize = "small")

## Figure A7: Conditional Average Treatment Effects on Moderates ----
### Create an empty data frame to store the estimates for later visualization
df.CATE1 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE1) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE1[1:3, 1] <- 1
df.CATE1[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (7-point scale DV)
### All respondents (including inattentive respondents)
H1a.full <- glht(mod3.moder, linfct = c("treatment1 + treatment1:moderate = 0"))
df.CATE1[1, 2] <- confint(H1a.full)$confint[ , "Estimate"]
df.CATE1[1, 3] <- confint(H1a.full)$confint[ , "lwr"]
df.CATE1[1, 4] <- confint(H1a.full)$confint[ , "upr"]
df.CATE1[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3.alt1 <- lm_robust(support ~ treatment1 + treatment2 + moderate +
                         treatment1 * moderate + treatment2 * moderate +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                       data = subset(df, attention == 1 | attention == 2))
H1a.alt1 <- glht(mod3.alt1, linfct = c("treatment1 + treatment1:moderate = 0"))
df.CATE1[2, 2] <- confint(H1a.alt1)$confint[ , "Estimate"]
df.CATE1[2, 3] <- confint(H1a.alt1)$confint[ , "lwr"]
df.CATE1[2, 4] <- confint(H1a.alt1)$confint[ , "upr"]
df.CATE1[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3.alt2 <- lm_robust(support ~ treatment1 + treatment2 + moderate +
                         treatment1 * moderate + treatment2 * moderate +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                       data = subset(df, attention == 2))
H1a.alt2 <- glht(mod3.alt2, linfct = c("treatment1 + treatment1:moderate = 0"))
df.CATE1[3, 2] <- confint(H1a.alt2)$confint[ , "Estimate"]
df.CATE1[3, 3] <- confint(H1a.alt2)$confint[ , "lwr"]
df.CATE1[3, 4] <- confint(H1a.alt2)$confint[ , "upr"]
df.CATE1[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (7-point scale DV)
### All respondents (including inattentive respondents)
H2a.full <- glht(mod3.moder, linfct = c("treatment2 + treatment2:moderate = 0"))
df.CATE1[4, 2] <- confint(H2a.full)$confint[ , "Estimate"]
df.CATE1[4, 3] <- confint(H2a.full)$confint[ , "lwr"]
df.CATE1[4, 4] <- confint(H2a.full)$confint[ , "upr"]
df.CATE1[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.alt1 <- glht(mod3.alt1, linfct = c("treatment2 + treatment2:moderate = 0"))
df.CATE1[5, 2] <- confint(H2a.alt1)$confint[ , "Estimate"]
df.CATE1[5, 3] <- confint(H2a.alt1)$confint[ , "lwr"]
df.CATE1[5, 4] <- confint(H2a.alt1)$confint[ , "upr"]
df.CATE1[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.alt2 <- glht(mod3.alt2, linfct = c("treatment2 + treatment2:moderate = 0"))
df.CATE1[6, 2] <- confint(H2a.alt2)$confint[ , "Estimate"]
df.CATE1[6, 3] <- confint(H2a.alt2)$confint[ , "lwr"]
df.CATE1[6, 4] <- confint(H2a.alt2)$confint[ , "upr"]
df.CATE1[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE1$group <- factor(df.CATE1$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE_moder1 <- 
  ggplot(data = df.CATE1, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Moderates") +
  ggtitle("Panel A: Seven-Point UBI Support") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.title = element_blank(),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-1.15, 1.15)) +
  guides(color = guide_legend(nrow = 1, byrow = T))
fig_CATE_moder1

### Create an empty data frame to store the estimates for later visualization
df.CATE2 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE2) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE2[1:3, 1] <- 1
df.CATE2[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (binary DV)
### All respondents (including inattentive respondents)
mod3.binary0 <- lm_robust(percent.sup ~ treatment1 + treatment2 + moderate +
                            treatment1 * moderate + treatment2 * moderate +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = df)
H1a.binary0 <- glht(mod3.binary0, linfct = c("treatment1 + treatment1:moderate = 0"))
df.CATE2[1, 2] <- confint(H1a.binary0)$confint[ , "Estimate"]
df.CATE2[1, 3] <- confint(H1a.binary0)$confint[ , "lwr"]
df.CATE2[1, 4] <- confint(H1a.binary0)$confint[ , "upr"]
df.CATE2[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3.binary1 <- lm_robust(percent.sup ~ treatment1 + treatment2 + moderate +
                            treatment1 * moderate + treatment2 * moderate +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = subset(df, attention == 1 | attention == 2))
H1a.binary1 <- glht(mod3.binary1, linfct = c("treatment1 + treatment1:moderate = 0"))
df.CATE2[2, 2] <- confint(H1a.binary1)$confint[ , "Estimate"]
df.CATE2[2, 3] <- confint(H1a.binary1)$confint[ , "lwr"]
df.CATE2[2, 4] <- confint(H1a.binary1)$confint[ , "upr"]
df.CATE2[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3.binary2 <- lm_robust(percent.sup ~ treatment1 + treatment2 + moderate +
                            treatment1 * moderate + treatment2 * moderate +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = subset(df, attention == 2))
H1a.binary2 <- glht(mod3.binary2, linfct = c("treatment1 + treatment1:moderate = 0"))
df.CATE2[3, 2] <- confint(H1a.binary2)$confint[ , "Estimate"]
df.CATE2[3, 3] <- confint(H1a.binary2)$confint[ , "lwr"]
df.CATE2[3, 4] <- confint(H1a.binary2)$confint[ , "upr"]
df.CATE2[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (binary DV)
### All respondents (including inattentive respondents)
H2a.binary0 <- glht(mod3.binary0, linfct = c("treatment2 + treatment2:moderate = 0"))
df.CATE2[4, 2] <- confint(H2a.binary0)$confint[ , "Estimate"]
df.CATE2[4, 3] <- confint(H2a.binary0)$confint[ , "lwr"]
df.CATE2[4, 4] <- confint(H2a.binary0)$confint[ , "upr"]
df.CATE2[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.binary1 <- glht(mod3.binary1, linfct = c("treatment2 + treatment2:moderate = 0"))
df.CATE2[5, 2] <- confint(H2a.binary1)$confint[ , "Estimate"]
df.CATE2[5, 3] <- confint(H2a.binary1)$confint[ , "lwr"]
df.CATE2[5, 4] <- confint(H2a.binary1)$confint[ , "upr"]
df.CATE2[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.binary2 <- glht(mod3.binary2, linfct = c("treatment2 + treatment2:moderate = 0"))
df.CATE2[6, 2] <- confint(H2a.binary2)$confint[ , "Estimate"]
df.CATE2[6, 3] <- confint(H2a.binary2)$confint[ , "lwr"]
df.CATE2[6, 4] <- confint(H2a.binary2)$confint[ , "upr"]
df.CATE2[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE2$group <- factor(df.CATE2$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE_moder2 <- 
  ggplot(data = df.CATE2, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Moderates") +
  ggtitle("Panel B: Percentage Support for UBI") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-.3, .3)) +
  guides(colour = guide_legend(nrow = 1, byrow = TRUE))
fig_CATE_moder2

### Combine the two plots (fig_CATE_self1 and fig_CATE_self2) and let them share a common legend
legend <- get_legend(fig_CATE_moder1)  # save the legend
fig_CATE_moderate <- 
  grid.arrange(arrangeGrob(fig_CATE_moder1 + theme(legend.position = "none"),
                           fig_CATE_moder2 + theme(legend.position = "none"),
                           nrow = 1),
               legend, nrow = 2, heights = c(10, 1))
ggsave(file = "fig_CATE_moderate.pdf", fig_CATE_moderate, width = 8, height = 4.5)

# Appendix G ----
## Figure A16: Heterogeneous Treatment Effects by Education Level ----
HTE_educ <- interflex(Y = "support", D = "group", X = "educ", 
                      Z = c("age", "female", "married", "educ", "income", "white",
                      "black", "employ", "union", "pid", "welfare.receive"), 
                      data = df, na.rm = T, estimator = "kernel")
HTE_educ_p1 <- 
  plot(HTE_educ, order = c("1", "2"), 
     subtitles = c("Equalizing Opportunity", "Limiting Government"),
     theme.bw = T, cex.sub = 1.2, cex.lab = .8, cex.axis = .8,
     xlab = "Moderator: Education Level (0 = Less than High School; 7 = Doctorate Degree)", 
     ylab = "Marginal Framing Effect on UBI Support")
ggsave(file = "fig_HTE_educ1.pdf", HTE_educ_p1, width = 8, height = 4.5)
HTE_educ_p2 <- 
  predict(HTE_educ, order = c("0", "1", "2"), 
     subtitles = c("Control Group", "Equalizing Opportunity", "Limiting Government"),
     pool = F, color = c("Salmon", "#7FC97F", "#BEAED4"), 
     ylim = c(0, 6), theme.bw = T, cex.sub = 1.2,
     xlab = "Moderator: Education Level (0 = Less than High School; 7 = Doctorate Degree)", 
     ylab = "Expected Level of UBI Support")
ggsave(file = "fig_HTE_educ2.pdf", HTE_educ_p2, width = 8, height = 4.5)

## Figure A17: Heterogeneous Treatment Effects by Income Level ----
HTE_income <- interflex(Y = "support", D = "group", X = "income", 
                        Z = c("age", "female", "married", "educ", "income", "white",
                              "black", "employ", "union", "pid", "welfare.receive"), 
                        data = df, na.rm = T, estimator = "kernel")
HTE_income_p1 <- 
  plot(HTE_income, order = c("1", "2"), 
       subtitles = c("Equalizing Opportunity", "Limiting Government"),
       xlim = c(0, 11.8), theme.bw = T, cex.sub = 1.2, cex.lab = .8, cex.axis = .6,
       xlab = "Moderator: Income Level (0 = Lowest; 12 = Highest)", 
       ylab = "Marginal Framing Effect on UBI Support")
ggsave(file = "fig_HTE_income1.pdf", HTE_income_p1, width = 8, height = 4.5)
HTE_income_p2 <- 
  predict(HTE_income, order = c("0", "1", "2"), 
          subtitles = c("Control Group", "Equalizing Opportunity", "Limiting Government"),
          pool = F, color = c("Salmon", "#7FC97F", "#BEAED4"), 
          xlim = c(0, 11.8), ylim = c(0, 6), theme.bw = T, cex.sub = 1.2, cex.axis = .8,
          xlab = "Moderator: Income Level (0 = Lowest; 12 = Highest)", 
          ylab = "Expected Level of UBI Support")
ggsave(file = "fig_HTE_income2.pdf", HTE_income_p2, width = 8, height = 4.5)